This script contains the example analysis presented during the Tumors and the immune system 2022 - High dimensional spatial profiling of tumor microenvironment course. It focuses on highly multiplexed image analysis after image processing and is a shortened version of the IMC data analysis book.

Getting started

To follow the tutorial, please clone the repository:

git clone https://github.com/nilseling/MultiplexedImageAnalysis2022.git

or download the zipped version by clicking on Code and then Download ZIP.

The MultiplexedImageAnalysis.Rmd file contains runnable code to follow the tutorial.

Install required packages

Please install R and RStudio.

To follow the data analysis examples in this tutorial, you will need to install the following R packages from within R:

install.packages("BiocManager")
BiocManager::install(c("cytomapper", "imcRtools", "openxlsx", "stringr", 
                       "dittoSeq", "tidyverse", "bluster", "patchwork",
                       "viridis", "scater", "scuttle", "S4Vectors", "igraph"))

Please also install FIJI and QuPATH.

Download example data

To highlight the basic steps of multiplexed image analysis, we provide example data that were acquired as part of the Integrated iMMUnoprofiling of large adaptive CANcer patient cohorts projects (immucan.eu). The raw data of 4 patients can be accessed online at zenodo.org/record/5949116. Image processing and segmentation has been performed using steinbock and we will only need to access the sample/patient information here:

download.file("https://zenodo.org/record/5949116/files/sample_metadata.xlsx", 
         destfile = "data/sample_metadata.xlsx")

Processed multiplexed imaging data

The IMC raw data was processed using the steinbock framework. Image processing included file type conversion, cell segmentation and feature extraction.

steinbock output

The needed output of the steinbock framework includes the single-cell mean intensity files, the single-cell morphological features and spatial locations, spatial object graphs in form of edge lists indication cells in close proximity, hot pixel filtered multi-channel images, segmentation masks, image metadata and channel metadata. All these files will be downloaded here for easy access. The commands which were used to process the data can be found in data/steinbock/steinbock.sh.

# download intensities
url <- "https://zenodo.org/record/6043600/files/intensities.zip"
destfile <- "data/steinbock/intensities.zip"
download.file(url, destfile)
unzip(destfile, exdir="data/steinbock", overwrite=TRUE)
unlink(destfile)

# download regionprops
url <- "https://zenodo.org/record/6043600/files/regionprops.zip"
destfile <- "data/steinbock/regionprops.zip"
download.file(url, destfile)
unzip(destfile, exdir="data/steinbock", overwrite=TRUE)
unlink(destfile)


# download neighbors
url <- "https://zenodo.org/record/6043600/files/neighbors.zip"
destfile <- "data/steinbock/neighbors.zip"
download.file(url, destfile)
unzip(destfile, exdir="data/steinbock", overwrite=TRUE)
unlink(destfile)

# download images
url <- "https://zenodo.org/record/6043600/files/img.zip"
destfile <- "data/steinbock/img.zip"
download.file(url, destfile)
unzip(destfile, exdir="data/steinbock", overwrite=TRUE)
unlink(destfile)

# download masks
url <- "https://zenodo.org/record/6043600/files/masks_deepcell.zip"
destfile <- "data/steinbock/masks_deepcell.zip"
download.file(url, destfile)
unzip(destfile, exdir="data/steinbock", overwrite=TRUE)
unlink(destfile)

# download individual files
download.file("https://zenodo.org/record/6043600/files/panel.csv", 
              "data/steinbock/panel.csv")
download.file("https://zenodo.org/record/6043600/files/images.csv", 
              "data/steinbock/images.csv")
download.file("https://zenodo.org/record/6043600/files/steinbock.sh", 
              "data/steinbock/steinbock.sh")

Visualization of multiplexed images

The first part of the tutorial will present FIJI and QuPATH for interactive multiplexed image visualization. Please make sure you have this software installed.

Handling image and single-cell data in R

This section describes how to read in single-cell data and images into R after image processing and segmentation (see Section @ref(processing)).

To highlight examples for IMC data analysis, we provide already processed data at https://zenodo.org/record/6043600. This data has already been downloaded in Section @ref(download-data) and can be accessed in the folder data.

We use the imcRtools package to read in single-cell data extracted using the steinbock framework or the IMC Segmentation Pipeline. Both image processing approaches also generate multi-channel images and segmentation masks that can be read into R using the cytomapper package.

library(imcRtools)
## Warning: package 'S4Vectors' was built under R version 4.1.3
library(cytomapper)

Reading in single-cell data

For single-cell data analysis in R the SingleCellExperiment [@Amezquita2019] data container is commonly used within the Bioconductor framework. It allows standardized access to (i) expression data, (ii) cellular metadata (e.g. cell type), (iii) feature metadata (e.g. marker name) and (iv) experiment-wide metadata. For an in-depth introduction to the SingleCellExperiment container, please refer to the SingleCellExperiment class.

The SpatialExperiment class [@Righelli2021] is an extension of the SingleCellExperiment class. It was developed to store spatial data in addition to single-cell data and an extended introduction is accessible here.

To read in single-cell data generated by the steinbock framework or the IMC Segmentation Pipeline, the imcRtools package provides the read_steinbock and read_cpout function, respectively. By default, the data is read into a SpatialExperiment object; however, data can be read in as a SingleCellExperiment object by setting return_as = "sce". All functions presented in this book are applicable to both data containers.

The steinbock framework was used to process and segment IMC raw data available here. In Section @ref(download-data), the processed data was downloaded from here and moved to the data/steinbock folder.

The read_steinbock function provided by imcRtools can now be used to read in steinbock generated data. For more information, please refer to ?read_steinbock.

spe <- read_steinbock("data/steinbock/")
spe
## class: SpatialExperiment 
## dim: 40 43289 
## metadata(0):
## assays(1): counts
## rownames(40): MPO HistoneH3 ... DNA1 DNA2
## rowData names(11): channel name ... Final.Concentration...Dilution
##   uL.to.add
## colnames: NULL
## colData names(8): sample_id ObjectNumber ... width_px height_px
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialData names(0) :
## spatialCoords names(2) : Pos_X Pos_Y
## imgData names(1): sample_id

By default, single-cell data is read in as SpatialExperiment object. The summarized pixel intensities per channel and cell (here mean intensity) are stored in the counts slot. Columns represent cells and rows represent channels.

counts(spe)[1:5,1:5]
##                 [,1]      [,2]      [,3]      [,4]      [,5]
## MPO        0.5000000 0.3750000 0.3435238 0.3407180 0.3348934
## HistoneH3 11.6130212 3.1264606 4.8283998 6.7566772 6.1166809
## SMA        1.8240418 0.2831075 0.9293683 0.9744561 5.7186161
## CD16       3.2254507 0.8346569 1.0149211 6.9618460 0.5949180
## CD38       0.8098145 0.4273930 0.4318003 1.1435687 0.3998025

Metadata associated to individual cells are stored in the colData slot. After initial image processing, these metadata include the numeric identifier (ObjectNumber), the area, and morphological features of each cell. In addition, sample_id stores the image name from which each cell was extracted and the width and height of the corresponding images are stored.

head(colData(spe))
## DataFrame with 6 rows and 8 columns
##      sample_id ObjectNumber      area major_axis_length minor_axis_length
##    <character>    <numeric> <numeric>         <numeric>         <numeric>
## 1 Patient1_001            1        22          12.64911           2.00000
## 2 Patient1_001            2         8           4.47214           2.00000
## 3 Patient1_001            3        28          10.67014           3.52327
## 4 Patient1_001            4        23           8.87999           3.21096
## 5 Patient1_001            5        27          10.73261           3.13466
## 6 Patient1_001            6        22           9.82805           2.88133
##   eccentricity  width_px height_px
##      <numeric> <numeric> <numeric>
## 1     0.987421       600       600
## 2     0.894427       600       600
## 3     0.943911       600       600
## 4     0.932335       600       600
## 5     0.956397       600       600
## 6     0.956059       600       600

The main difference between the SpatialExperiment and the SingleCellExperiment data container in the current setting is the way spatial locations of all cells are stored. For the SingleCellExperiment container, the locations are stored in the colData slot while the SpatialExperiment container stores them in the spatialCoords slot:

head(spatialCoords(spe))
##      Pos_X     Pos_Y
## 1 517.0000 0.5000000
## 2 231.5000 0.5000000
## 3 271.4643 1.0000000
## 4 287.3478 0.9565217
## 5 410.6296 0.8888889
## 6 468.7273 0.7727273

The spatial object graphs generated by steinbock (see Section @ref(feature-extraction) are read into a colPair slot of the SpatialExperiment (or SingleCellExperiment) object. Cell-cell interactions (cells in close spatial proximity) are represented as “edge list” (stored as SelfHits object). Here, the left side represents the column indices of the “from” cells and the right side represents the column indices of the “to” cells. For visualization of the spatial object graphs, please refer to Section @ref(spatial-viz).

colPair(spe, "neighborhood")
## SelfHits object with 202720 hits and 0 metadata columns:
##                 from        to
##            <integer> <integer>
##        [1]         1        41
##        [2]         1        69
##        [3]         2        18
##        [4]         2        58
##        [5]         2        85
##        ...       ...       ...
##   [202716]     43288     43249
##   [202717]     43288     43268
##   [202718]     43288     43269
##   [202719]     43289     43248
##   [202720]     43289     43266
##   -------
##   nnode: 43289

Finally, metadata regarding the channels are stored in the rowData slot. This information is extracted from the panel.csv file. Channels are ordered by isotope mass and therefore match the channel order of the multi-channel images (see Section @ref(read-images)).

head(rowData(spe))
## DataFrame with 6 rows and 11 columns
##               channel        name      keep   ilastik  deepcell Tube.Number
##           <character> <character> <numeric> <numeric> <numeric>   <numeric>
## MPO               Y89         MPO         1        NA        NA        2101
## HistoneH3       In113   HistoneH3         1         1         1        2113
## SMA             In115         SMA         1        NA        NA        1914
## CD16            Pr141        CD16         1        NA        NA        2079
## CD38            Nd142        CD38         1        NA        NA        2095
## HLADR           Nd143       HLADR         1        NA        NA        2087
##                        Target Antibody.Clone Stock.Concentration
##                   <character>    <character>           <numeric>
## MPO       Myeloperoxidase MPO Polyclonal MPO                 500
## HistoneH3          Histone H3           D1H2                 500
## SMA                       SMA            1A4                 500
## CD16                     CD16       EPR16784                 500
## CD38                     CD38        EPR4106                 500
## HLADR                  HLA-DR        TAL 1B5                 500
##           Final.Concentration...Dilution   uL.to.add
##                              <character> <character>
## MPO                              4 ug/mL         0.8
## HistoneH3                        1 ug/mL         0.2
## SMA                           0.25 ug/mL        0.05
## CD16                             5 ug/mL           1
## CD38                           2.5 ug/mL         0.5
## HLADR                            1 ug/mL         0.2

After reading in the single-cell data, few further processing steps need to be taken.

Add additional metadata

We can set the colnames of the object to generate unique identifiers per cell:

colnames(spe) <- paste0(spe$sample_id, "_", spe$ObjectNumber)

It is also often the case that sample-specific metadata are available externally. For the current data, we need to link the cancer type (also referred to as “Indication”) to each sample. This metadata is available as external excel file:

library(openxlsx)
library(stringr)
meta <- read.xlsx("data/sample_metadata.xlsx")
spe$patient_id <- as.vector(str_extract_all(spe$sample_id, "Patient[1-4]", simplify = TRUE))
spe$ROI <- as.vector(str_extract_all(spe$sample_id, "00[1-8]", simplify = TRUE))
spe$indication <- meta$Indication[match(spe$patient_id, meta$Sample.ID)]
unique(spe$indication)
## [1] "SCCHN" "BCC"   "NSCLC" "CRC"

The selected patients were diagnosed with different cancer types:

  • SCCHN - head and neck cancer
  • BCC - breast cancer
  • NSCLC - lung cancer
  • CRC - colorectal cancer

Transform counts

The distribution of expression counts across cells is often observed to be skewed towards the right side meaning lots of cells display low counts and few cells have high counts. To avoid analysis biases from these high-expressing cells, the expression counts are commonly transformed or clipped.

Here, we perform counts transformation using an inverse hyperbolic sine function. This transformation is commonly applied to flow cytometry data. The cofactor here defines the expression range on which no scaling is performed. While the cofactor for CyTOF data is often set to 5, IMC data usually display much lower counts. We therefore apply a cofactor of 1.

However, other transformations such as log(counts(spe) + 0.01) should be tested when analysing IMC data.

library(dittoSeq)
dittoRidgePlot(spe, var = "CD3", group.by = "patient_id", assay = "counts") +
    ggtitle("CD3 - before transformation")

assay(spe, "exprs") <- asinh(counts(spe)/1)
dittoRidgePlot(spe, var = "CD3", group.by = "patient_id", assay = "exprs") +
    ggtitle("CD3 - after transformation")

Define interesting channels

For downstream analysis such as visualization, dimensionality reduction and clustering, only a subset of markers should be used. As convenience, we can store an additional entry in the rowData slot that specifies the markers of interest. Here, we deselect the nuclear markers and keep all other biological targets.

rowData(spe)$use_channel <- !grepl("DNA|Histone", rownames(spe))

Reading in images

The cytomapper package allows multi-channel image handling and visualization within the Bioconductor framework. The most common data format for multi-channel images or segmentation masks is the TIFF file format, which is used by steinbock to save images.

Here, we will read in multi-channel images and segmentation masks into a CytoImageList data container. It allows storing multiple multi-channel images and requires matched channels across all images within the object.

The loadImages function is used to read in processed multi-channel images and their corresponding segmentation masks as generated by steinbock. The multi-channel images are saved as 32-bit images while the segmentation masks are saved as 16-bit images. To correctly scale pixel values of the segmentation masks when reading them in set as.is = TRUE.

images <- loadImages("data/steinbock/img/")
## All files in the provided location will be read in.
masks <- loadImages("data/steinbock/masks_deepcell/", as.is = TRUE)
## All files in the provided location will be read in.

In the case of multi-channel images, it is beneficial to set the channelNames for easy visualization. Using the steinbock framework, the channel order of the single-cell data matches the channel order of the multi-channel images. However, it is recommended to make sure that the channel order is known.

channelNames(images) <- rownames(spe)
images
## CytoImageList containing 14 image(s)
## names(14): Patient1_001 Patient1_002 Patient1_003 Patient2_001 Patient2_002 Patient2_003 Patient2_004 Patient3_001 Patient3_002 Patient3_003 Patient4_005 Patient4_006 Patient4_007 Patient4_008 
## Each image contains 40 channel(s)
## channelNames(40): MPO HistoneH3 SMA CD16 CD38 HLADR CD27 CD15 CD45RA CD163 B2M CD20 CD68 Ido1 CD3 LAG3 / LAG33 CD11c PD1 PDGFRb CD7 GrzB PDL1 TCF7 CD45RO FOXP3 ICOS CD8a CarbonicAnhydrase CD33 Ki67 VISTA CD40 CD4 CD14 Ecad CD303 CD206 cleavedPARP DNA1 DNA2

For visualization shown in Section @ref(image-visualization) we will need to add additional metadata to the elementMetadata slot of the CytoImageList objects. This slot is easily accessible using the mcols function.

Here, we will save the matched sample_id, patient_id and indication information within the elementMetadata slot of the multi-channel images and segmentation masks objects. It is crucial that the order of the images in both CytoImageList objects is the same.

all.equal(names(images), names(masks))
## [1] TRUE
patient_id <- str_extract_all(names(images), "Patient[1-4]", simplify = TRUE)
indication <- meta$Indication[match(patient_id, meta$Sample.ID)] 
mcols(images) <- mcols(masks) <- DataFrame(sample_id = names(images),
                                           patient_id = patient_id,
                                           indication = indication)

Single-cell analysis

The single-cell data contained in the SpatialExperiment object can now be analysed to extract biological knowledge. As part of the first exploratory analysis steps, we want to visualize the cell density per mm\(^2\).

library(tidyverse)

colData(spe) %>%
    as.data.frame() %>%
    group_by(sample_id) %>%
    summarize(cell_count = n(),
           no_pixels = mean(width_px) * mean(height_px)) %>%
    mutate(cells_per_mm2 = cell_count/(no_pixels/1000000)) %>%
    ggplot() +
        geom_point(aes(sample_id, cells_per_mm2)) + 
        theme_minimal(base_size = 15)  + 
        theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
        ylab("Cells per mm2") + xlab("")

A convenient way to visualize the protein expression in individual cells is using heatmaps. In the following code chunk, we will randomly select 2000 cells and visualize their expression.

library(dittoSeq)
library(viridis)
cur_cells <- sample(seq_len(ncol(spe)), 2000)
dittoHeatmap(spe[,cur_cells], genes = rownames(spe)[rowData(spe)$use_channel],
             assay = "exprs", cluster_cols = TRUE, scale = "none",
             heatmap.colors = viridis(100), annot.by = "indication")

An alternative way of visualizing high-dimensional data is to first perform dimwnsionality reduction followed by scatter plot visualization of the cells in 2 dimensions. A common approach to perform non-linear dimensionality reduction is Uniform Manifold Approximation and Projection (UMAP). The scater Bioconductor package can be used to perform this dimensionality reduction.

library(scater)
set.seed(220404)
spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs") 

After dimensionality reduction, the low-dimensional embeddings are stored in the reducedDim slot.

head(reducedDim(spe, "UMAP"))
##                     [,1]      [,2]
## Patient1_001_1 -4.087332 -3.196866
## Patient1_001_2 -5.046047 -3.088559
## Patient1_001_3 -4.762129 -3.672599
## Patient1_001_4 -3.937544 -2.925449
## Patient1_001_5 -1.389223  3.402379
## Patient1_001_6 -4.377478 -3.541066

We can now visualize cells in low dimensions, color them by patient id and visualize protein expression.

library(patchwork)

# visualize patient id 
p1 <- dittoDimPlot(spe, var = "patient_id", reduction.use = "UMAP", size = 0.2) +
    ggtitle("Patient ID on UMAP")

# visualize region of interest id
p2 <- dittoDimPlot(spe, var = "ROI", reduction.use = "UMAP", size = 0.2) +
    ggtitle("ROI ID on UMAP")

p1 + p2

# visualize marker expression
p1 <- dittoDimPlot(spe, var = "Ecad", reduction.use = "UMAP", 
                   assay = "exprs", size = 0.2) +
    scale_color_viridis(name = "Ecad") +
    ggtitle("E-Cadherin expression on UMAP")
p2 <- dittoDimPlot(spe, var = "CD45RO", reduction.use = "UMAP", 
                   assay = "exprs", size = 0.2) +
    scale_color_viridis(name = "CD45RO") +
    ggtitle("CD45RO expression on UMAP")

p1 + p2

A common step in single-cell data analysis is the identification of cellular phenotypes. This is often done by clustering cells based on their similarity in protein expression. For clustering, we will use the bluster package and specifically using a graph-based clustering approach.

library(bluster)
mat <- t(assay(spe, "exprs")[rowData(spe)$use_channel,])

cluster <- clusterRows(mat, NNGraphParam(k=30, cluster.fun = "louvain"))
length(unique(cluster))
## [1] 14
spe$cluster <- as.factor(cluster)

Using this approach, we detect 13 clusters. To annotate clusters, we will need to visualize the protein expression per cluster. For this, we will visualize a heatmap of the 2000 samples cells and observe single-cell cluster expression.

Based on the expression, the individual clusters can then be annotated.

dittoHeatmap(spe[,cur_cells], genes = rownames(spe)[rowData(spe)$use_channel],
             assay = "exprs", cluster_cols = FALSE, annot.by = "cluster",
             scaled.to.max = TRUE, heatmap.colors.max.scaled = viridis(100))

cluster_annot <- as.numeric(cluster)
cluster_annot[cluster_annot == 1] <- "Tumor_proliferating"
cluster_annot[cluster_annot == 2] <- "Stroma"
cluster_annot[cluster_annot == 3] <- "Tumor_proliferating"
cluster_annot[cluster_annot == 4] <- "Stroma"
cluster_annot[cluster_annot == 5] <- "CD38"
cluster_annot[cluster_annot == 6] <- "Macrophages"
cluster_annot[cluster_annot == 7] <- "Tumor_hypoxic"
cluster_annot[cluster_annot == 8] <- "helper_Tcell"
cluster_annot[cluster_annot == 9] <- "cyto_Tcell"
cluster_annot[cluster_annot == 10] <- "Stroma"
cluster_annot[cluster_annot == 11] <- "Tumor"
cluster_annot[cluster_annot == 12] <- "Tumor_proliferating"
cluster_annot[cluster_annot == 13] <- "B_cells"

spe$cluster_annot <- as.factor(cluster_annot)

We can now also visualize the location and interactions between the identified cell types in the tissue. For this, we select randomly 4 images.

set.seed(1234)
cur_img_ids <- sample(unique(spe$sample_id), 4)

plotSpatial(spe[,spe$sample_id %in% cur_img_ids], img_id = "sample_id", 
            node_color_by = "cluster_annot", draw_edges = TRUE, colPairName = "neighborhood",
            nodes_first = FALSE)

This concludes the single-cell analysis and visulization and we will now move on to the image-level visualization.

Image visualization

In the next section, we will visualize some of the images of the dataset. For this, we will use the cytomapper package. We will consistently visualize the four images randomly sampled above.

Pixel level visualization

The cytomapper package can be used to visualize pixel-level information by merging channels and therefore generating image composites that visualize more than one single protein.

In the first example, we will vizualice E-cadherin (tumor cells), Ki67 (proliferation marker), CD3 (T cells) and CD20 (B cells) and a nuclear marker.

plotPixels(images[cur_img_ids],
           colour_by = c("Ecad", "Ki67", "CD3", "CD20", "DNA1"),
           colour = list(Ecad = c("black", "brown"),
                         Ki67 = c("black", "yellow"),
                         CD3 = c("black", "green"),
                         CD20 = c("black", "red"),
                         DNA1 = c("black", "blue")),
           bcg = list(
               Ecad = c(0, 10, 1),
               Ki67 = c(0, 10, 1),
               CD3 = c(0, 10, 1),
               CD20 = c(0, 15, 1),
               DNA1 = c(0, 5, 1)
           ),
           legend = list(colour_by.title.cex = 0.7,
                         colour_by.labels.cex = 0.7))

On these images we can now outline selected cell types. This approach is useful to visualize if the corrected cells were phenotyped. In the following example, T cells are outlined.

plotPixels(images[cur_img_ids],
           mask = masks[cur_img_ids],
           object = spe[,spe$cluster_annot == "cyto_Tcell"],
           img_id = "sample_id",
           cell_id = "ObjectNumber",
           outline_by = "cluster_annot",
           missing_colour = "white",
           thick = TRUE,
           colour_by = c("Ecad", "Ki67", "CD8a", "CD20", "DNA1"),
           colour = list(Ecad = c("black", "brown"),
                         Ki67 = c("black", "yellow"),
                         CD8a = c("black", "green"),
                         CD20 = c("black", "red"),
                         DNA1 = c("black", "blue")),
           bcg = list(
               Ecad = c(0, 10, 1),
               Ki67 = c(0, 10, 1),
               CD8a = c(0, 10, 1),
               CD20 = c(0, 15, 1),
               DNA1 = c(0, 5, 1)
           ),
           legend = list(colour_by.title.cex = 0.7,
                         colour_by.labels.cex = 0.7))

### Cell level visualization

Next to the pixel-level visualization, the cytomapper package provides the plotCells function to visualize cell-level metadata on segmentation masks.

Here, we can visualize the defined cell phenotypes on segmentation masks.

plotCells(mask = masks[cur_img_ids],
          object = spe,
          img_id = "sample_id",
          cell_id = "ObjectNumber",
          colour_by = "cluster_annot")

This sort of visualization is also useful when selecting images that contain certain tissue structures:

plotCells(mask = masks,
          object = spe,
          img_id = "sample_id",
          cell_id = "ObjectNumber",
          colour_by = "cluster_annot")

Software used

The following chunk lists all software packages used in this tutorial.

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] bluster_1.4.0               patchwork_1.1.1            
##  [3] scater_1.22.0               scuttle_1.4.0              
##  [5] viridis_0.6.2               viridisLite_0.4.0          
##  [7] forcats_0.5.1               dplyr_1.0.8                
##  [9] purrr_0.3.4                 readr_2.1.2                
## [11] tidyr_1.2.0                 tibble_3.1.6               
## [13] tidyverse_1.3.1             dittoSeq_1.7.0             
## [15] ggplot2_3.3.5               stringr_1.4.0              
## [17] openxlsx_4.2.5              cytomapper_1.7.1           
## [19] EBImage_4.36.0              imcRtools_1.1.8            
## [21] SpatialExperiment_1.4.0     SingleCellExperiment_1.16.0
## [23] SummarizedExperiment_1.24.0 Biobase_2.54.0             
## [25] GenomicRanges_1.46.1        GenomeInfoDb_1.30.1        
## [27] IRanges_2.28.0              S4Vectors_0.32.4           
## [29] BiocGenerics_0.40.0         MatrixGenerics_1.6.0       
## [31] matrixStats_0.61.0         
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2                shinydashboard_0.7.2     
##   [3] R.utils_2.11.0            tidyselect_1.1.2         
##   [5] htmlwidgets_1.5.4         grid_4.1.2               
##   [7] BiocParallel_1.28.3       DropletUtils_1.14.2      
##   [9] ScaledMatrix_1.2.0        munsell_0.5.0            
##  [11] codetools_0.2-18          units_0.8-0              
##  [13] DT_0.22                   withr_2.5.0              
##  [15] colorspace_2.0-3          highr_0.9                
##  [17] knitr_1.38                rstudioapi_0.13          
##  [19] labeling_0.4.2            GenomeInfoDbData_1.2.7   
##  [21] polyclip_1.10-0           bit64_4.0.5              
##  [23] farver_2.1.0              pheatmap_1.0.12          
##  [25] rhdf5_2.38.1              vctrs_0.4.0              
##  [27] generics_0.1.2            xfun_0.30                
##  [29] R6_2.5.1                  ggbeeswarm_0.6.0         
##  [31] graphlayouts_0.8.0        rsvd_1.0.5               
##  [33] locfit_1.5-9.5            concaveman_1.1.0         
##  [35] bitops_1.0-7              rhdf5filters_1.6.0       
##  [37] RTriangle_1.6-0.10        DelayedArray_0.20.0      
##  [39] assertthat_0.2.1          promises_1.2.0.1         
##  [41] scales_1.1.1              vroom_1.5.7              
##  [43] ggraph_2.0.5              beeswarm_0.4.0           
##  [45] gtable_0.3.0              beachmat_2.10.0          
##  [47] tidygraph_1.2.0           rlang_1.0.2              
##  [49] systemfonts_1.0.4         broom_0.7.12             
##  [51] yaml_2.3.5                abind_1.4-5              
##  [53] modelr_0.1.8              backports_1.4.1          
##  [55] httpuv_1.6.5              tools_4.1.2              
##  [57] ellipsis_0.3.2            raster_3.5-15            
##  [59] jquerylib_0.1.4           RColorBrewer_1.1-3       
##  [61] proxy_0.4-26              ggridges_0.5.3           
##  [63] Rcpp_1.0.8.3              plyr_1.8.7               
##  [65] sparseMatrixStats_1.6.0   zlibbioc_1.40.0          
##  [67] classInt_0.4-3            RCurl_1.98-1.6           
##  [69] cowplot_1.1.1             cluster_2.1.2            
##  [71] haven_2.4.3               ggrepel_0.9.1            
##  [73] fs_1.5.2                  magrittr_2.0.3           
##  [75] RSpectra_0.16-0           data.table_1.14.2        
##  [77] magick_2.7.3              reprex_2.0.1             
##  [79] hms_1.1.1                 mime_0.12                
##  [81] fftwtools_0.9-11          evaluate_0.15            
##  [83] xtable_1.8-4              jpeg_0.1-9               
##  [85] readxl_1.4.0              gridExtra_2.3            
##  [87] compiler_4.1.2            KernSmooth_2.23-20       
##  [89] crayon_1.5.1              R.oo_1.24.0              
##  [91] htmltools_0.5.2           later_1.3.0              
##  [93] tzdb_0.3.0                tiff_0.1-11              
##  [95] lubridate_1.8.0           DBI_1.1.2                
##  [97] tweenr_1.0.2              dbplyr_2.1.1             
##  [99] MASS_7.3-54               sf_1.0-7                 
## [101] Matrix_1.4-1              cli_3.2.0                
## [103] R.methodsS3_1.8.1         parallel_4.1.2           
## [105] igraph_1.3.0              pkgconfig_2.0.3          
## [107] sp_1.4-6                  terra_1.5-21             
## [109] xml2_1.3.3                svglite_2.1.0            
## [111] vipor_0.4.5               bslib_0.3.1              
## [113] dqrng_0.3.0               XVector_0.34.0           
## [115] rvest_1.0.2               digest_0.6.29            
## [117] RcppAnnoy_0.0.19          rmarkdown_2.13           
## [119] cellranger_1.1.0          uwot_0.1.11              
## [121] edgeR_3.36.0              DelayedMatrixStats_1.16.0
## [123] shiny_1.7.1               rjson_0.2.21             
## [125] lifecycle_1.0.1           jsonlite_1.8.0           
## [127] Rhdf5lib_1.16.0           BiocNeighbors_1.12.0     
## [129] limma_3.50.1              fansi_1.0.3              
## [131] pillar_1.7.0              lattice_0.20-45          
## [133] fastmap_1.1.0             httr_1.4.2               
## [135] glue_1.6.2                zip_2.2.0                
## [137] png_0.1-7                 svgPanZoom_0.3.4         
## [139] bit_4.0.4                 ggforce_0.3.3            
## [141] class_7.3-19              stringi_1.7.6            
## [143] sass_0.4.1                HDF5Array_1.22.1         
## [145] nnls_1.4                  BiocSingular_1.10.0      
## [147] irlba_2.3.5               e1071_1.7-9